On the anomalous dynamics of capillary rise in porous media 
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The anomalous dynamics of capillary rise in a porous medium discovered experimentally more 
than a decade ago (Delker et al., Phys. Rev. Lett. 76 (1996) 2902) is described. The developed 
theory is based on considering the principal modes of motion of the menisci that collectively form 
the wetting front on the Darcy scale. These modes, which include (i) dynamic wetting mode, (ii) 
threshold mode and (iii) interface de-pinning process, are incorporated into the boundary conditions 
for the bulk equations formulated in the regular framework of continuum mechanics of porous media, 
thus allowing one to consider a general case of three-dimensional flows. The developed theory makes 
it possible to describe all regimes observed in the experiment, with the time spanning more than 
four orders of magnitude, and highlights the dominant physical mechanisms at different stages of 
the process. 
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I. INTRODUCTION 

The propagation of the liquid-fluid interfaces through 
porous media is central to a wide range of natural 
phenomena and industrial applications, with the lat- 
ter including enhanced oil recovery, hydrogeology, fuel 
cells, carbon dioxide sequestration to mention but a 
few. This topic remains the subject of intensive re- 
search, both experimental and theoretical, comprehen- 
sively reviewed in a number of articles over the past forty 
years 0, H, [H, Hoi HH- The main aspects of research 
have been the rate of propagation of the wetting front 
[H [H, [H UpL the wetting front's roughening (jj, HzJ 
and stability [20|, as well as the related problems of the 
formation and dynamics of the pockets of the displaced 
phase (bubbles, ganglia) left behind the front [il.l20l.l28l]. 
The first of these aspects is of particular importance as it 
ultimately determines the main macroscopic characteris- 
tics of the process in many applications. 

As has been discovered experimentally more than a 
decade ago by Delker and his co-workers ||, besides 
the common situation where a wetting front propagates 
through a porous medium broadly in accordance with 
Washburn's model [3fj| . which balances the driving force 
due to the (presumed constant) capillary pressure of a 
meniscus and viscous resistance as in the Poiseuille-type 
flow, for some media, such as porous matrices made of 
packed spherical beads, the initial Washburn-type imbi- 
bition is followed by a completely different and in many 
ways 'anomalous' regime. A representative set of data 
taken from Q is shown in Fig. [TJ Similar results have 
been reported later by Lago and Araujo [TtJ • The essence 



of the discovered effect is that, if the height h of the cap- 
illary rise (measured from some initial level to remove 
from consideration the entrance effects) is plotted against 
time t on the log- log scale (Fig. [JJ, one can immediately 
see two distinct regions. Roughly speaking, for about 
two minutes the liquid climbs 2/3 of its eventual (max- 
imum) height h m ax m the Washburn-like regime after 
which it takes many hours for it to advances across the 
remaining 1/3 of h max , with the wetting front moving 
in small-amplitude jumps on the pore scale [17[- The 
log-log plot of this second regime shows a clear concave- 
convex sequence which indicates that the dynamics there 
is more complex than what one would expect from some 
power-law fit and the accompanying arguments. Another 
intriguing feature of the phenomenon is that h max is de- 
termined by the balance of capillarity and gravity, i.e. the 
factors that, together with viscous resistance, determine 
the dynamics of the Washburn regime, although what 
looked like the Washburn regime has been abandoned af- 
ter a couple of minutes from the onset of the capillary rise 
and for hours the process is distinctly non-Washburnian. 

The experimental data has been discussed qualitatively 
in terms of interface pinning and random capillary forces 
[1, [TtJ , but on the quantitative level the only outcome is 
that a simple equation 



dh/dt = v (F/F T - 1) 



/3 



(1) 



expressing the rate of the capillary rise as a function of 
a driving force F and a threshold value Ft leads to an 
"anomalously large" exponent ft Q, so that h deviates 
from the experimental data for small times and unphys- 
ically diverges as time goes to infinity. The dashed line 
in Fig. [TJ corresponds to 



h = H c -(H c -hi)[l+A(t-t 1 )] 



V(i- 



(2) 
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that has been deduced from ([TJ and used in ||; the val- 
ues of the constants H c , hi, A, ti and are given in 0]. 
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Although the fitting curve ([2]) is able to describe only a 
finite time span, the general ideas of the interface pin- 
ning and random forces that might de-pin the interface 
and allow it to move further seem fruitful, and a question 
that arises naturally is how to embed them into the regu- 
lar framework of continuum mechanics of porous media, 
as opposed to just using ad hoc semi-empirical equations 
for the wetting front evolution in a one-dimensional flow. 
Below, we address this question on the basis of an ear- 
lier developed approach to the modelling of the wetting 
front dynamics in porous media based on considering dif- 
ferent modes of motion that menisci go through on the 
pore scale and the corresponding technique of conjugate 
problems [26]. Then, we will discuss how the experi- 
mental phenomenon in question is seen through other 
modelling approaches. 

II. MACROSCOPIC (DARCY-SCALE) 
DESCRIPTION 

In the continuum framework, the wetting front dQ\ is a 
moving boundary which, together with other boundaries 
dfl2, confines a domain f2 where the Darcy equation, 

u = -{K/n)V(p + pgz), (re SI), (3) 

and the continuity equation, V-u = 0, for the average ve- 
locity u and pressure p operate. Then, the wetting front 
evolution is part of the solution of a properly formulated 
problem for these bulk equations. The Darcy equation 
([3]) is written in the form already accounting for gravity 
with p being the density of the liquid, g the gravitational 
acceleration, z the coordinate directed against gravity, n 
the permeability of the porous matrix and \x the liquid's 
viscosity; the coordinates are represented in terms of the 
position vector r; hereafter the pressure is measured with 
respect to the (presumed constant) pressure in the dis- 
placed gas. 

An appropriate starting point for the modelling is the 
recently developed approach [2(| that gives boundary 
conditions for Laplace's equation for p, 

V 2 p = 0, (re Q), (4) 

which follows from the Darcy and continuity equations 
above and can be used to replace the latter. The key 
idea of this approach is to consider the modes of mo- 
tion which the menisci that collectively form the wet- 
ting front undergo as the wetting front propagates. The 
simplest model formulated in the framework of this ap- 
proach accounts for the two main modes: (i) the wetting 
mode, where, on the pore scale, the contact line moves 
forward, essentially in the Washburn regime but account- 
ing for a dynamic, i.e. velocity-dependent, contact angle 
9d, and (ii) the threshold mode, where the contact line 
gets pinned and the meniscus bends as the pressure on 
it increases until the contact angle reaches some thresh- 
old value when the meniscus can go back into the 








FIG. 1. Time-dependence of the imbibition height in ID 
capillary rise. Circles: experimental data from [8| for beads 
253 (jm in diameter; dashed line: fit considered in [s|; solid 
line: present theory. 

wetting mode with the contact line moving again. This 
increase of pressure on the meniscus as the contact line 
gets pinned is similar to what one would have on a piston 
sucking a liquid into a pipe if the motion of the piston is 
blocked. For the porous medium, the maximum possible 
pressure on the meniscus in the threshold mode, p\dfin 
is the solution of a conjugate problem [26j 

V 2 p = 0, (refi); n V(p + pgz)\ dni = 0, (5) 

with the boundary condition for p on dQ 2 being the same 
as for p; n is the outward normal to dSli. Thus, the idea 
of the interface pinning is already in the model, used 
in formulating the boundary conditions on the wetting 
front that we will recapitulate below, and in this work we 
consider how the model can be generalized to incorporate 
the idea of random forces that could lead to de-pinning 
of the interface and describe the observed features of the 
phenomenon mentioned earlier. 

On the moving wetting front the kinematic and dy- 
namic boundary conditions for (j4]) have the form 

^ + u-V/ = 0, (6) 



P = A lPl +A 2 p 2 , (reffii), (7) 

where for a one-dimensional capillary rise /(r, t) = z — 
h(t), px, p2 are the averaged pressures and Ai, An. are 
the spatio-temporally averaged fractions of the unit area 
of the free surface corresponding to the two modes of 
motion {A\ + A 2 = 1). For the wetting mode one has 

pi = —2oco&9d/a, (8) 

where a is the liquid-gas surface tension, a is the effective 
radius of the capillary and the dependence of the dynamic 
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contact angle 9d on the meniscus speed u\ is given by 

u%_ = / (! + (1 -pf e )cosg fl )(cosg 8 -cose d ) 2 \ 1/2 

Uci V 4(coBfl,+B)(cQsfl«j + B) J ' 

(9) 

where £> = (1 — p\ e ) l (\ + p s le uo(8d)), 9 S is the static 
contact angle, 



sin 9 d - d cos 9 d /7pg(l + 4a/3)\ 

?Wd) = -t— 7, 7. tt, U c i = 

sin 6 d cos 9d-0d \ t/3 J 



1/2 



is the characteristic speed associated with the param- 
eters that the 'additional' physics of wetting brings in 
to resolve the well-known 'moving contact-line problem' 
[TIT l24l . [25j . Pq, p\ e , a, (3, 7, t are material constants 
characterizing the contacting media whose values can be 
found elsewhere @, [HJ • 

In the threshold mode, the contact line gets pinned and 
the meniscus, experiencing an increase in pressure on it, 
bends so that the contact angle varies from 8 d , i.e. the 
value with which the meniscus arrives at the threshold 
mode, to which is the value at which the meniscus 
'breaks through' the threshold and the process goes back 
to the wetting mode. From the dynamics of this type of 
motion, under the assumption that the meniscus retains 
the shape of a spherical cap with its radius varying as 
the meniscus bends, one has that [261 



u 2 = aJ(8 d )T~ 



where 



(10) 



P2 = P~^ h^+cos0 d \J(0 d ), (11) 



(12) 



2a 



Pa 
2a~ 



J(9 d ) 
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Plan!, 
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fQ 7T 




- tan 
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V2 ~ 4 





tan 



[f] a = f(b) — /(a), and the time T that the meniscus 
spends in the threshold mode is given by 



t = t 2 e d ,— \ = 



Pa 
2a~ 



Pa\ _ 



h cos 9 C 

2a 



I[9 d ,— 
1 ' 2a 



(13) 



'2a) Jg d (1 + sin 9) 2 (Pa /{2a) + cos0)' 

For a one-dimensional capillary rise, as follows from ([5]). 
the stagnation pressure p\dn 1 is given simply by pld^ = 
Po — pgh(t), where po is the prescribed pressure at z = 0. 

Finally, the velocity of the wetting front as a whole, 
U n = n ■ ujao^ is given by an expression 



which is similar to 0. For the menisci intermittently 
going through the wetting and threshold modes as the 
wetting front propagates, the coefficients A\ and A2 can 
be viewed as reflecting the fraction of the time spent in 
each mode as the meniscus travels over the length of av- 
eraging that introduces the Darcy scale (or, equivalently, 
the fraction of the interfacial area corresponding to each 
mode of motion over the time interval of averaging that 
introduces the Darcy time scale), i.e. the spatio-temporal 
averages we mentioned earlier. They are given by [26J 



Ay = 



SiU 2 



-4, 



S 2 Ul 



S2U1 + S\U 2 



S 2 Ui + SiU 2 



(15) 



where 



«i( 







1. 

sio, 



9 d - 9, > 
9 d - 9, < 



S2 



(16) 



AyUi + A 2 U 2 , 



(14) 



and sio (< 1) is a characteristic of the porous matrix. 
Then, as should be expected, the slowest (controlling) 
mode of motion makes a greater contribution to the av- 
erage pressure and velocity at the wetting front, and if the 
velocity ui corresponding to the ith mode reaches zero, 
one will have A.- L — 1 and the pressure at the wetting 
front, that is now at rest, will become equal to pi. Thus, 
the wetting front will stop propagating in two cases: (a) 
u\ = and hence 6 d — 9 Sl which means that the menis- 
cus has reached its equilibrium state corresponding to the 
maximum imbibition height h max , and (b) u 2 — so that 
the wetting front still has a capacity to propagate further 
but the contact line became pinned (threshold mode) and 
the pressure that mounts on the meniscus in this case, 
even when it reaches its maximum possible value p\orin 
is insufficient to push the meniscus through. Mathemat- 
ically, in the last case we have that p\dn ± , which goes 
down as h increases, becomes equal to = — 2crcos#*/a 
and hence is unable to make the contact angle greater 
than 0*, which would allow the meniscus to resume its 
motion in the wetting mode: in this case I(9d, Pa/ (2a)) 
and hence the time T go to infinity. The height corre- 
sponding to this last case depends only on 6**, and the 
meniscus reaches it in a finite time 1261. 



III. SUBCRITICAL INTERFACE DE-PINNING 

It is convenient to introduce the 'stagnation' contact 
angle 9 corresponding to the stagnation pressure p|ar2i by 
9 = arccos(— p\gQ 1 a/(2aj). Then, if 9 > 9 S , the wetting 
front has the potential to propagate but if at the same 
time 9 < 9* the stagnation pressure p\dn t is unable to 
push the meniscus through the threshold mode. In a 
sense, 9* — 9 can be viewed as a quantitative measure 
of the potential barrier that has to be overcome to get 
the meniscus back into the wetting mode when plasii is 
'subcritical', i.e. less than p*, and hence unable to push 
the meniscus through the threshold. 

Importantly, since, until the wetting front reaches its 
equilibrium position at the maximum height, at every 
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moment individual menisci are not in the same mode of 
motion (and, for the threshold motion, not even in the 
same stage of it), the Darcy-scale pressures we are consid- 
ering, including the stagnation pressure £>|ani> represent 
average values, whereas on the pore scale one also has 
pressure fluctuations due to mutual influences of menisci. 
These fluctuations are unimportant when the stagnation 
pressure p\dn 1 is capable of pushing the meniscus through 
the threshold mode. However, as p\dn 1 goes down to p*, 
the time T needed to overcome the threshold increases, so 
that, when it becomes large enough, it is the fluctuations 
that increasingly become the mechanism of de-pinning, 
and when 9 < 6* it is only the random fluctuations that 
can de-pin the menisci. 

The simplest way of accounting for the subcritical 
de-pinning due to random fluctuations as this mecha- 
nism takes over from the 'regular' de-pinning due to the 
stagnation pressure is to assume that, once T 2 becomes 
greater than a certain value T+, it is these random fac- 
tors that will de-pin the interface and determine the time 
it stays in the threshold mode. For the 'regular' stagna- 
tion pressure one has that T 2 — > 00 as p\dQ t —> P* 
and p\dn t is no longer able to push the meniscus through 
when p\dn 1 < P*- The probability of the random fac- 
tors de-pinning the interface should be expected to go 
down as p* — p\dQ t (or equivalently (9* — 9) increases. 
The simplest way of generalizing the model to incorpo- 
rate the above scenario mathematically is to replace (fT"3]) 
and (|12p respectively with 

T=< U (Od, » if ^(planj < T+ ^ ^ 

[T+ + k(9+ - 9) 2 , if T 2 (p\ dni ) > T+ 

P _ Jplenx, if T 2 {p\aa 1 ) <T+ , , 

\P+, XT 2 (p\ dni )>T + ' W 

where p + is determined by T 2 (p + ) = 2+, and 9 + = 
arccos(— p+a/ (2a)). Since T 2 rises steeply only when 9 is 
close to 9* , in practice one has that 9+ « 9* and p+ 

Now, we have a closed model, which, unlike ad hoc 
formulae for one-dimensional propagation of the wetting 
front, is applicable for a general case of three-dimensional 
flows. In order to describe a particular flow involving a 
moving wetting front, one has to solve Laplace's equation 
((4]) for p in fi whose boundary d£li evolves according 
to (|6]), where u is given by ((SJ, subject to the dynamic 
condition 0, where equations ([8])-(jTTj). (fl ^ -(flg ]| close 
the formulation and p\dOi is the solution of the conjugate 
problem ([5]); the boundary conditions on dfl 2 for both p 
and p are the same and, together with the initial shape 
of f2, they specify a particular problem. 

In the case of a one-dimensional capillary rise Laplace's 
equations for p and p give that these are linear functions 
of z, and equations (J3]), Q, © yield 

dh k { pq — p(h, t) \ 
dt fx \ h ^ J ' 



which together with the algebraic equations ([7]), where 
p\ 9Ql =p(h,t), ©-d}, d, where u n = dh/dt, fity- 
(|18p . with plant — Po — pgh a s the solution of the con- 
jugate problem, form a closed system for h, p(h,t), pi, 
9d, Ui, u 2 , p 2 , A\, A 2 , s\ and s 2 . The results of com- 
paring the numerical solution corresponding to the ex- 
perimental flow conditions of [8] with the data are shown 
in Fig. [TJ As one can see, the solid curve representing 
the computed solution describes the data very well over 
the whole time period spanning more than four orders 
of magnitude. The theoretical curve also levels off as 
t — > 00, indicating that the capillary rise does eventu- 
ally come to a halt. Comparison of the theory with all 
four sets of experimental data from Q is shown in Fig. [2] 
(The dashed line in this figure is used to indicate that for 
the beads' diameter of 510 /im, strictly speaking, the the- 
ory is used outside its limits of applicability as the whole 
advancement of the wetting front is less than 40 beads' di- 
ameters, so that it is difficult to talk about the separation 
of scales required for the continuum mechanics approach 
to work. Indeed, for this approach to be applicable, there 
should exist an intermediate scale much larger than the 
pore size and at the same time much smaller than the 
macroscopic length scale on which the flow is described. 
In the case of 510 /im beads, 'much large' and 'much 
smaller' would mean 6 times larger or smaller, which is 
clearly not sufficient to ensure acceptable accuracy.) 

It is noteworthy that, although the initial regime where 
the curve in Fig. [T] rises steeply looks Washburn-likc, 
it actually involves both the wetting and the thresh- 
old modes of motion. As in (2(|, the presence of the 
threshold-overcoming motion becomes pronounced only 
when h climbs close to = (po —p*)/(pg), i.e. when 
p becomes close to p* or, in other words, 9 close to 0*. 
For the results presented in Fig. [1] 9* = 67°, 9 S = 0°, 
sio = 0.7, fiU c i/(Kpg) = 10 2 , p{ e = 0.6, T+/T = 3, 
k/T — 4 x 10 3 , where T = 2a/j,/((pg) 2 aK), and it is 
9 S , 9*, T + /Tq and k/T^ that are most important. Since 
in the experiment, as described in 0, the bottom of the 
test section was located approximately 4 cm above the 
bottom of the porous sample, we have to start the calcu- 
lations from the bottom of the sample with pq = pgZi n i, 
Zini = 4 cm and, to compare theory with the data, mea- 
sure the time from the moment the wetting front reaches 
Zi n i. To describe the data in Fig. [2 only a variation of 
9* and K = k/T is required: 6**i 8 o = 64°, 6>* 2 53 = 67°, 
0*359 = 78°; K 359 = 3K 253 = 12K 18Q ., K 180 = 10 3 . 

It is worth mentioning that, although in the model we 
represent capillary effects in the pores using the capillary 
pressure and viscous resistance corresponding to an 'ef- 
fective' circular cross-section, whereas in the experiment 
with the porous matrix made of spherical beads neither 
the 'chambers' nor the 'throats' of the porous medium 
had circular cross-sections, no adjustment of the results 
was required: we used the radius of the bead as a in the 
model with no subsequent calibration of the time and 
length scales. This shows that a 'representative' way 
of modelling the porous medium, as opposed to much 
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FIG. 2. All four experimental sets from [8j for beads' diame- 
ters 180 ^m (V), 253 /im (o), 359 (im (a), and 510 (□). 
The corresponding theoretical curves are given as solid lines 
for the first 3 sets of data and as a dashed line for the last. 

more difficult way of incorporating into a model the exact 
porous structure determined via elaborate and expensive 
experiments, allows one to incorporate all the main fea- 
tures of the process on the pore scale, including the actual 
physics of dynamic wetting, and obtain good results for 
the flow on the Darcy scale. 



IV. DISCUSSION 

It is instructive to look at the obtained theoretical re- 
sult and the experiment it addresses from the viewpoint 
of the different modelling approaches used to describe 
two-phase flows in porous media. These approaches 
broadly fall into two general classes, 'representative' and 
'simulative'. They are not antagonistic as, in theory, if 
the same pore-scale physics and the same characteristics 
of the porous medium are accounted for in the models 
formulated in the framework of each of these approaches, 
then the results produced by these models are expected 
to converge and describe the same macroscopic behaviour 
of the wetting front. 

The present model has been developed in the frame- 
work of the 'representative' approach, where the equa- 
tions and boundary conditions arc formulated on the 
Darcy scale, i.e. the scale implying that the continuum 
limit has already been taken, and the properties of the 
porous matrix are 'represented' in terms of the perme- 
ability coefficient (or tensor, if the porous medium is 
anisotropic), effective size of the pores (or the corre- 
sponding distribution), effective threshold angles, etc. 
Importantly, since the pore-scale wetting process is mod- 
elled realistically, with the velocity (as well as material) 
dependence of the dynamic contact angle, the model cap- 
tures naturally that the wetting front has the capacity to 
move forward when the contact angle is greater than the 
static value 9 S , i.e. when the interface has not reached its 



maximum height determined by the balance of capillary 
and gravity forces. In other words, as the meniscus gets 
de-pinned, the fact that the contact angle differs from its 
equilibrium value 6 S and that the dynamic contact angle 
is velocity-dependent act as a mechanism that moves the 
interface until the dynamic contact angle goes down to 
S and the interface reaches its maximum height. 

The subcritical de-pinning mechanism that comes into 
action when p|ani < V* is formulated implicitly, in terms 
of the 'potential barrier' 6 + — 6 and the 'waiting time' 
T + + k{6 + — 9) 2 required for the random fluctuations to 
overcome it. Both of these characteristics are Darcy-scale 
parameters. A natural way to develop the model further 
would be to remove the direct link between the 'poten- 
tial barrier' and the 'waiting time' and instead explic- 
itly introduce the field of pressure fluctuations depend- 
ing on the flow rate and properties of the porous matrix. 
Then, this pressure fluctuation field becomes an addition 
to p\dQi a $ a breakthrough factor. For this explicitly in- 
troduced mechanism, the results obtained in the present 
work would serve as a guideline, indicating one of the 
outcomes that this mechanism should produce. 

The implicit mechanism of the subcritical interface de- 
pinning on the Darcy scale that we have introduced can 
be viewed as a macroscopic manifestation of the dynam- 
ics of avalanches [l(| ■ Avalanches qualitatively stem from 
the same physics as the one considered here and, in a 
sense, can be regarded as a medium-scale phenomenon, 
i.e. between the pore scale and the Darcy scale. The 
idea of linking the Darcy-scale subcritical interface de- 
pinning and the dynamics of avalanches agrees with the 
fact that the avalanche-type events have been observed 
experimentally by Lago and Araujo [l7j in the anomalous 
regime of the wetting front propagation. The avalanche 
behaviour of the wetting fronts is currently being investi- 
gated in terms of scaling laws [22| , and the present theory, 
with its explicit accounting for different pressures corre- 
sponding to different modes of motion and their spatio- 
temporal weights, offers a macroscopic framework for the 
mathematical description of these medium-scale events. 
Then, the subcritical interface de-pinning in the anoma- 
lous regime of the wetting front propagation could be 
viewed as the macroscopic outcome of a succession of 
avalanches with decreasing probabilities. 

The 'simulative' approach to the modelling of two- 
phase flows in porous media is based on replacing the ac- 
tual porous matrix with a regular network of 'chambers' 
and capillary 'throats' connecting them pi. IH FP3, [TBL [l8j . 
In order to partially compensate the anisotropy inher- 
ent in this approach where, as is the case in most works, 
the chambers are placed at the nodes of a regular lat- 
tice, the sizes of both chambers and throats are made 
random following certain prescribed distributions. The 
macroscopic (Darcy-scale) characteristics of the process 
are obtained as a result of the appropriate averaging. The 
simulative approach has the appeal of what looks like a 
numerical experiment offering a transparent link between 
the pore-scale and the Darcy-scale dynamics, but, unlike 
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the case of molecular-dynamics simulations with regard 
to macroscopic fluid mechanics, this appeal is moderated 
by a number of factors, notably the fact that the actual 
dynamics on the pore scale is not computed. Instead, 
it is essentially represented in terms of a Washburn-type 
dynamics, thus by-passing the 'moving contact-line prob- 
le m' II II] and the associated physics of dynamic wetting 
@,H1[- It is also important to note that the rigid geomet- 
ric structure of the network simulating the actual porous 
matrix imposes unavoidable limitations on the macro- 
scopic transport properties of the porous medium that 
the network is purports to simulate. By contrast, the 
representative approach can introduce any tensorial and 
topological characteristics of the porous medium that the 
corresponding model requires. Having in mind the above 
shortcomings of the representative approach, it is inter- 
esting to look at what it can produce with regard to the 
anomalous regime of the wetting front propagation. 

Bijeljic et al. Q performed the network modelling of 
the capillary rise and compared their results to the trun- 
cated set of data taken from Lago and Araujo [l7|, with 
the original experiment by Delker et al. [8| mentioned 
but not used for comparison with the simulations. The 
simulations agree well with the experimental data for the 
times t < 4 x 10 4 s after which the simulated height levels 
off. However, when one takes the full set of experimental 
data reported by Lago and Araujo, i.e. for the times up 
to t = 2 x 10 5 s, one can immediately see that the wet- 
ting front continues to climb. It is also worth noting that 
the simulated height- vs-time curve in the log-log coordi- 
nates is convex, with the slope monotonically decreasing, 
whereas, as one can see in Fig. [1] the experimental data 
show a distinct concave-convex sequence, i.e. after an ini- 
tial decrease as the anomalous regime is entered the slope 
picks up again until, finally, the data start to level off, 
asymptotically approaching the maximum height. The 
same trend was observed earlier by Diggins et al. M 
whose data are given in Fig. 12 of Lago and Araujo [17j . 
This figure in [l7| also clearly shows that, although the 
time-dependence of the height of the capillary rise follow- 
ing from the Washburn-type interplay of capillarity, vis- 
cosity and gravity can describe the 'regular' regime and 
be fitted to the initial stage of the 'anomalous' regime, 
it is nowhere near a satisfactory description of the latter 
once the full set of data is considered, as the height-vs- 
time curve in the logarithmic coordinates picks up again 
and climbs much higher than the above fitting predicts. 

From the viewpoint of the theory developed in the 
present work, the main deficiency of the currently imple- 
mented network models is that they essentially deal only 
with one — wetting — mode of motion of the menisci. 
Then, setting aside the minor (in comparison with the ef- 
fects considered here) variations introduced by randomiz- 



ing the size distributions, these models broadly reproduce 
the Washburn-type dynamics of the wetting front. As a 
result, once gravity starts to balance capillarity as the 
driving force, the wetting front slows down and comes 
to a halt. Essentially, the fitting of the simulated curve 
to the experimental data for the very beginning of the 
anomalous regime is produced by adjusting the maximum 
height of the capillary rise, whereas, as we pointed out in 
the introduction, the intriguing feature of the anomalous 
regime is precisely the fact that it lies in between the 
'normal' Washburnian regime of the imbibition and the 
maximum height that is also Washburnian. 

The capillary network approach can be modified in a 
relatively straightforward way to account for the thresh- 
old mode of motion. The main element in this modifica- 
tion should be equipping the 'chambers' with threshold 
characteristics, such as the pressure required to overcome 
the threshold, which, besides material properties, can de- 
pend on the number of menisci reaching the same cham- 
ber. The implementation of the subcritical de-pinning 
is more challenging as this would require accounting for 
fluctuations of pressure experienced by the liquid, i.e. 
replacing the Washburn-type models of the flow in the 
'throats' by an essentially unsteady motion that takes 
into account the unsteady processes in the neighboring 
chambers and throats. An intermediate check for such a 
model could be its ability to produce avalanches as the 
medium-scale phenomena that on the Darcy scale result 
in the anomalous regime of imbibition. 

V. CONCLUSION 

The developed theory shows that the new approach 
to the modelling of the propagation of wetting fronts in 
porous media based on considering specific modes of mo- 
tion that the menisci of the pore scale undergo as the 
front propagates allows one to incorporate critical phe- 
nomena and adequately describe experimental data for 
the anomalous regime of imbibition. Accounting for the 
random pore-scale forces macroscopically, in terms of the 
'potential barriers' and the corresponding times required 
for the random forces to overcome these barriers, allowed 
the simplest model formulated in the framework of the 
new approach to describe the whole experimental curve, 
from the Washburn regime to the (also Washburn) max- 
imum imbibition height with the anomalous regime in 
between. The proposed theory could be used as a guide 
for the porous network modelling and the study of the 
anomalous imbibition regime as the manifestation of the 
dynamics of avalanches. 

This publication was based on work supported in part 
by Award No. KUK-C1-013-04 , made by King Abdullah 
University of Science and Technology (KAUST). 
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